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Abstract 



A theoretical analysis is given of the equation of motion method, due to 
Alben et al., to compute the eigenvalue distribution (density of states) of 
very large matrices. The salient feature of this method is that for matrices of 
the kind encountered in quantum physics the memory and CPU requirements 
of this method scale linearly with the dimension of the matrix. We derive a 
rigorous estimate of the statistical error, supporting earlier observations that 
the computational efficiency of this approach increases with matrix size. We 
use this method and an imaginary-time version of it to compute the energy 
and the specific heat of three different, exactly solvable, spin-1/2 models and 
compare with the exact results to study the dependence of the statistical 
errors on sample and matrix size. 



1 



Typeset using REVT^ 



I. INTRODUCTION 



The calculation of the distribution of eigenvalues of very large matrices is a central 
problem in quantum physics. This distribution determines the thermodynamic properties 
of the system (see below). It is directly related to the single-particle density of states (DOS) 
or Green's function. In a one-particle (e.g., one-electron) description knowledge of the DOS 
suffices to compute the transport properties 

The most direct method to compute the DOS, i.e. all the eigenvalues, is to diagonalize 
the matrix H representing the Hamiltonian of the system. This approach has two obvious 
limitations: The number of operations increases as the third power of the dimension D 
of H and, perhaps most importantly, the amount of memory required by state-of-the-art 
algorithms grows as 0)0] • This scaling behavior limits the application of this approach 
to matrices of dimension D = (9(10000), which is too small for many problems of interest. 
What is needed are methods that scale linearly with D. 

There has been considerable interest in developing "fast" (i.e. 0{D)) algorithms to 
compute the DOS and other similar quantities. One such algorithm and an application 
of it to electron motion in disordered alloy models was given by Alben et al. In this 
approach the DOS is obtained by solving the time-dependent Schrodinger equation (TDSE) 
of a particle moving on a lattice, followed by a Fourier transform of the retarded Green's 
function Q. Using the unconditionally stable split-step Fast Fourier Transform (FFT) 
method to solve the TDSE, it was shown that the eigenvalue spectrum of a particle moving 
in continuum space can be computed in the same manner Fast algorithms of this 
kind proved useful to study various aspects of localization of waves §-§1 and other one- 
particle problems PJTO|]. Apphcation of these ideas to quantum many-body systems triggered 
further development of flexible and efficient methods to solve the TDSE. Based on Suzuki's 
product formula approach, an unconditionally stable algorithm was developed and used to 
compute the time-evolution of two-dimensional S=l/2 Heisenberg-like models |Tl|]. Results 
for the DOS of matrices of dimension D ^ 1000000 where reported |Tl|]. A potentially 
interesting feature of these fast algorithms is that they may run very efficiently on a quantum 
computer |T^,[TB[ . 



A common feature of these fast algorithms is that they solve the TDSE for a sample 
of randomly chosen initial states. The efficiency of this approach as a whole relies on 
the hypothesis (suggested by the central limit theorem) that satisfactory accuracy can be 
achieved by using a small sample of initial states. Experience not only shows that this 
hypothesis is correct, it strongly suggests that for a fixed sample size the statistical error on 
physical quantities such as the energy and specific heat decreases with the dimension D of 
the Hilbert space \T2\. 



In view of the general applicability of these fast algorithms to a wide variety of quantum 
problems it seems warranted to analyze in detail their properties and the peculiar D depen- 
dence in particular. In Sections and |TT| we recapitulate the essence of the approach. We 
present a rigorous estimate for the mean square error (variance) on the trace of a matrix. In 
Section |IV| we describe the imaginary-time version of the method. The statistical analysis of 
the numerical data is discussed in Section 0. Section ^ describes the model systems that 
are used in our numerical experiments. The algorithm used to solve the TDSE is reviewed 
in Section [VI I|. In Section [VIII| we derive rigorous bounds on the accuracy with which all 
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eigenvalues can be determined and demonstrate that this accuracy decreases hnearly with 
the time over which the TDSE is solved. The results of our numerical calculations are 
presented in Section ^ and our conclusions are given in Section 



II. THEORY 

The trace of a matrix A acting on a dimensional Hilbert-space spanned by an or- 
thonormal set of states {|0n)} is given by 

Tr A= X:(<^n|A<^n). (1) 
n=l 

Note that according to (1) we have Tt 1 = D. If D is very large one might think of approx- 
imating Eq. (1) by sampling over a subset of K {K -C D) "important" basis vectors. The 
problem with this approach is that the notion "important" may be very model-dependent. 
Therefore it is better to sample in a different manner. We construct a random vector lip) by 
choosing D complex random numbers, Cn = fn + ig-n, with mean 0, for n = 1 . . . D, so 

|^) = Ec„|0„), (2) 

n=l 

and calculate 

{ij\Ai:) = J2 CcMA(Pn). (3) 

n,m=l 

If we now sample over S realizations of the random vectors {ip} and calculate the average, 
we obtain 

^ s s D 

^E(V^pI^V^p) = oE E CpCn,p{<Pm\A(Pn). (4) 
p=l p=l n,m=l 

Assuming that there is no correlation between the random numbers in different realizations 
and that the random numbers p and gn,p are drawn from an even and symmetric (both 
with respect to each variable) probabihty distribution (see Appendix for more details), 
we have 

1 ^ 

J™„ 9 ^ ^*m,P^n,p = E (|cp) (5) 
p=l 

where E {.) denotes the expectation value with respect to the probability distribution used 
to generate the c„,p's. In the r.h.s of (5) the subscripts of c„_p have been dropped to indicate 
that the expectation value does not depend on n or p. It follows immediately that 

jim 1 f](V;p|A^p) = e{\c\')TtA = E (|cp) E {4>n\A<Pn), (6) 

p=l n=l 
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showing that we can compute the trace of A by samphng over random states {ipp}, provided 
there is an efficient algorithm to calculate {ijjp\Aijjp) (see Section VII ). 



According to the central limit theorem, for a large but finite S we have 

^ E 4,,c,,p = E {\c\') + O (^-^) , (7) 



p 



meaning that the statistical error on the trace vanishes like 1 / V^, which is not surprising. 
What is surprising is that one can prove a much stronger result as follows. Let us first 
normalize the c„_p's so that, for all p, 



D 



E K,p\' = 1. (8) 



n=l 

This innocent looking step has far reaching consequences. First we note that the normaliza- 
tion renders the method exact in (the rather trivial) case that the matrix A is proportional 
to the unit matrix. The price we pay for this is that for fixed p, the Cn^p are now correlated 
but that does not cause problems (see Appendix Second it follows that i?(|cp) = 1/D. 
Obviously the error can be written as 

TTA-^j2i^,\Ai;p) = TTRA, (9) 
p=i 

where 

D ^ 

p=l 

is a traceless (due to Eq. (8)) Hermitian matrix of random numbers. We put X = Ti RA 
and compute i?(|Xp). The result for the general case can be found in Appendix 0. For 
a uniform distribution of the c„p's on the hyper-sphere defined by Yln=i kn.pl^ = 1 the 
expression simplifies considerably and we find 

E ( Ti RAf) = ^ ^, 11) 

an exact expression for the variance in terms of the sample size S, the dimension D of the 
matrix A and the (unknown) constants TrA'^'A and | Try4|. 
Invoking a generalization of Markov's inequality |I^ 

P(|X|2 > a) < ^^^^ ; Va>0, (12) 

where P{Q) denotes the probability for the statement Q to be true. We find that the 
probability that | Tri?y4p exceeds a fraction a of | TrAp is bounded by 

/|Tri?Ap \ 1 DTtA^A-\TtA\^ , , 

> a < ' ^ ; Va>0, (13) 



TrA|2 - - aS{D\\) |TrA| 
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or, in other words, the relative statistical error ca on the estimator of the trace of A is given 
by 



ca = 



DTiA^A- |TrA|2 
\i S{D + 1)\TtA\^ ' 



(14) 



if I Tr A| > 0. We see that = if A is proportional to a unit matrix. From (14) it follows 
that, in general, we may expect to vanish with the square root of SD. The prefactor is 
a measure for the relative spread of the eigenvalues of A and is obviously model dependent. 
The dependence of ca on S, D and the spectrum of A is corroborated by the numerical 
results presented below. 

It is also of interest to examine the effect of not normalizing the c^^p's. A calculation 
similar to the one that lead to the above results yields 



eA 



^''■^•■^ (15) 



\ ^|TrA|2- 

Clearly this bound is less sharp and does not vanish if A is proportional to a unit matrix. 

III. REAL-TIME METHOD 

The distribution of eigenvalues or density of states (DOS) of a quantum system is defined 

as 

D 

271 



n=l 



oo 



Vie) = Y,5{e-E^) = ^ j Tr e"^*^ dt, (16) 



oo 



where H is the Hamiltonian of the system and n runs over all the eigenvalues of H. The 
DOS contains all the physical information about the equilibrium properties of the system. 
For instance the partition function, the energy, and the heat capacity are given by 

("OO 

Z= / deV[e)e-^\ (17) 



E = - 



1 

- deeV{e)e-^% (18) 

Zj J — oo 

C = p'{^yyee'V{e)e~^^-E'y (19) 

respectively. Here (3 = l/fc^T and ks is Boltzmann's constant (we put ks = 1 and h = 1 
from now on). 

As explained above the trace in the integral (|16D can be estimated by sampling over 
random vectors. For the statistical error analysis discussed below it is convenient to define 
a DOS-per-sample by 



1 



oo 



d,{e) = / e^*^ (^p|e-**^^p)cit, (20) 
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where the subscript p labels the particular realization of the random state Itpp). The DOS 
is then given by 

V{e)=\im ^j:d,{e). (21) 

^ p=i 

Schematically the algorithm to compute dp{e) consists of the following steps: 

1. Generate a random state lippiO)), set t = 0. 

2. Copy this state to \ipp(t)). 

3. Calculate {■ipp{0)\tpp{t)) and store the result. 

4. Solve the TDSE for a small time step r, replacing \ipp(t)) by \ipp{t + r)) (see Section |V11| 
for model specific details). 

5. Repeat times from Step 3. 

6. Perform a Fourier transform on the tabulated result and store dp{e). 

In practice the Fourier transform in Eq. (0) is performed by the Fast Fourier Transform 
(FFT). We use a Gaussian window to account for the finite time tN used in the numerical 
time-integration of the TDSE. The number of time step determines the accuracy with 
which the eigenvalues can be computed. In Section VIII we prove that this systematic error 
in the eigenvalues vanishes as 1/tN. 

Since for any reasonable physical system (or finite matrix) the smallest eigenvalue Eq 
is finite, for all practical purposes dp{e) = for e < eo < -Eq- The value of eo is easily 
determined by examination of the bottom of spectrum. To compute Z, E, or C we simply 
replace the interval [— oo, +oo] by [eo,+C)o]. 



IV. IMAGINARY-TIME METHOD 

The real-time approach has the advantage that it yields information on all eigenvalues 
and can be used to compute both dynamic and static properties without suffering from 
numerical instabilities. However for the computation of the thermodynamic properties, the 
imaginary-time version is more efficient. We will use the imaginary-time method as an 
independent check on the results obtained by the real-time algorithm. 

Repeating the steps that lead to Eq. (17) we find 

Z = Tiexp{-l3H) 
1 ^ 

= hm - expi-/3H)^p), (22) 

S p=i 

with similar expressions for E and C. 
Furthermore we have 

ii^plH^e-^''^,) = {e-^"l^i,,\H^e-^^/'ij,), (23) 
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assuming H is Hermitian. Therefore we only need to propagate the random state for an 
imaginary time (3/2 instead of j3. Furthermore we do not need to perform an FFT. Disre- 
garding these minor differences, the algorithm is the same as in the real-time case with r 
replaced by —it. 



V. ERROR ANALYSIS 

Estimating the statistical error on the partition function Z is easy because it depends 
linearly on the trace of the (imaginary) time evolution operator. However the error on E 
and C depends on this trace in a more complicated manner and this fact has to be taken 
into account. 



First we define 



roo 

Zp= / dedp{e)e-^\ (24) 

roo 

hp = / dedp{e)ee-'^', (25) 

roo 

Wp= / dedp{e)e'^e->^', (26) 



for the real-time method and 



z^ = (V'p|e-^^Vp), (27) 

hp = {iPp\He-^''iPp), (28) 

Wp = {'iljp\H''e-f'"ijp), (29) 

for the imaginary-time method. 

For each value of P we generate the data {zp}, {hp}, and {wp}, for p = 1, . . . , 5*. For 
both cases we have 

Z = hm z, (30) 

5— +00 

E = lim -, (31) 

where x = J2p=i Xp- The standard deviations on z, h, and w are given by 

- /l^. (33) 
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where var(a;) = — denotes the variance on the data {xp}. However the sets of data 
{zp}, {hp} and {wp} are correlated since they are calculated from the same set {iV'p)}- 
These correlations in the data are accounted for by calculating the covariance matrix Mk^i 
{k,l = 1, . . . ,3) the elements of which are given by XkXi — Xk xj, where {xi}, {X2}, and {x^} 
are a shorthand for {zp}, {hp}, and {wp} respectively. The estimates for the errors in Z, E 
and C are given by 

= ^5z^ (36) 

6E^ - ^ ^ M ^'^ (37) 

S 1 ^ dxi^ dxi 

S -l^^^'^'^ dx^dxl' ^^^^ 

2 

where E = h/'z and and C = (3'^{w/'z — h /^^). 



VI. EXACTLY SOLVABLE SPIN 1/2 MODELS 

The most direct way to assess the validity of the approach described above is to carry out 
numerical experiments on exactly solvable models. In this paper we consider three different 
exactly solvable models, two spin-1/2 chains and a mean-field spin-1/2 model. The former 
have a complicated spectrum, the latter has a highly degenerate eigenvalue distribution. 
These spin models differ from those studied elsewhere [llll , p!2[| in that they belong to the 
class of integrable systems. 



A. Spin chains 

Open spin chains of L sites described by the Hamiltonian 



L-l 



(39) 



i=l 



i=l 



where cxf , af , and af denote the Pauli matrices and J, A and h are model parameters, can 
be solved exactly. They can be reduced to diagonal form by means of the Jordan- Wigner 
transformation DTSf. We have 



+ hL, 



where cf and Cj are spin-less fermion operators and 

= -J(l + A){6,j^i + 6i^ij) - 2h6ij, 



B. 



-J(l-A)(5,,_i-(5,_i,, 



(40) 



(41) 
(42) 
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are Lx L matrices. By further canonical transformation this Hamiltonian can be written as 



^ ^ 1\ 1 



k=l 



H = Y,Ak[nk--j+-TTA + hL, (43) 



where Uk is the number operator of state k and the A^'s are given by the solution of the 
eigenvalue equation 

{A-B){A + B)(Pk = Al<Pk. (44) 

In the general case this eigenvalue problem of the L x L Hermitian matrix {A — B){A + B) 
is most easily solved numerically. In the present paper we confine ourselves to two limiting 
cases: The XY model (A = 1) and the Ising model in a transverse field (A = 0). 

B. Mean field model 



The Hamiltonian of the mean-field model reads 

L L 

L 



H = -jY: ^.■ff^-hY.^t, (45) 

i>j=l i=l 



and can be rewritten as 



with 



i=l 

The single spin-L/2 Hamiltonian has eigenvalues 



H = -2^S-S-2hS' + -J, (46) 
L 2 



S=lj:^^- (47) 



with degeneracy 



Ei^m = -2J/(/ + l)/L - 2hm + (48) 



"''-= L/2 + / + l iL/2-/j- (^9) 



This rather trivial model serves as a test for the case of highly degenerate eigenvalues. 



VII. TIME EVOLUTION 

For the approach outlined in Sections |ITT| and ^ to be of practical use it is necessary 
that the matrix elements of the exponential of H can be calculated efficiently. The purpose 
of this section is to describe how this can be done. 
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The general form of the Hamiltonians of the models we study is 



H 



a 



(50) 



i,j=l Q=x,y,z 



i=l a=x,y,z 



where the first sum runs over all pairs P of spins, af {a = x, y, z) denotes the a-th component 
of the spin-1/2 operator representing the z-th spin. For both methods, we have to calculate 
the evolution of a random state, i.e. U{t)\'iIj) = exp{—iTH)\ip) or U{t)\iP) = exp{—TH)\ip) 
for the real and imaginary time method respectively. We will discuss the real-time case only, 
the imaginary-time problem can be solved in the same manner. 

Using the semi-group property U(ti)U(t2) = U(ti + we can write U{t) = U{t)"^ 
where t = mr. Then the main step is to replace U{t) by a symmetrized product-formula 
approximation . For the case at hand it is expedient to take 



f/(r) ^ f/(r) =e 



-iTH^/2 -iTHy/2 -irH^ 



-iTHy/2 -iTH^/2 
^ 1 



(51) 



where 



E 

i=\ 



a = x,y, z. 



(52) 



Other decompositions |TI]JT^ work equally well but are somewhat less efficient for the cases 
at hand. In the real-time approach U{t) is unitary and hence the method is unconditionally 
stable |]T^ (also the imaginary-time method can be made unconditionally stable). It can 
be shown that ||f/(T) — f/(T)|| < sr^ (s > a constant) |jT9|, implying that the algorithm 
is correct to second order in the time step r [T^. Usually it is not difficult to choose r so 
small that for all practical purposes the results obtained can be considered as being "exact" . 
Moreover, if necessary, U{t) can be used as a building block to construct higher-order 
algorithms pO|-p^. In Appendix ^ we will derive bounds on the error in the eigenvalues 
when they are calculated using a symmetric product formula. 

As basis states {|0n)} we take the direct product of the eigenvectors of the S- (i.e. spin-up 
and spin-down lU)). In this basis, e"*"^^^/^ changes the input state by altering the phase 
of each of the basis vectors. As is a sum of pair interactions it is trivial to rewrite this 
operation as a direct product of 4x4 diagonal matrices (containing the interaction-controlled 
phase shifts) and 4x4 unit matrices. Still working in the same representation, the action of 
^-iTHy/2 ^g^^ written in a similar manner but the matrices that contain the interaction- 
controlled phase-shift have to be replaced by non-diagonal matrices. Although this does 
not present a real problem it is more efficient and systematic to proceed as follows. Let us 
denote by X{Y) the rotation by 7r/2 of each spin about the x(?/)-axis. As 



^-^rHy/2 ^ X e-'^^'y/^ X = Xe-""''^''' x\ 



(53) 



it is clear that the action of e can be computed by applying to each spin, the inverse of 

X followed by an interaction-controlled phase-shift and X itself. The prime in (K^ indicates 



that Jf and 



in Hz have to be replaced by J\j and }v{ respectively. A similar procedure 



is used to compute the action of e 



We only have to replace X hj Y. 
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VIII. ACCURACY OF THE COMPUTED EIGENVALUES 



First we consider the problem of how to choose the number of time steps to obtain 
the DOS with acceptable accuracy. According to the Nyquist sampling theorem employing 
a sampling interval At = n / maxj \Ei\ is sufficient to cover the full range of eigenvalues. On 
the other hand the time step also determines the accuracy of the approximation U{t). Let 
us call the maximum value of r which gives satisfactory accuracy tq (for the imaginary-time 
method, this is the only parameter). For the examples treated here tq < At), implying that 
we have to use more steps to solve the TDSE than we actually use to compute the FFT. 
Eigenvalues that differ less than Ae = ix/NAt cannot be identified properly. However since 
Ae oc A^^^ we only have to extend the length of the calculation by a factor of two to increase 
the resolution by the same factor. 

At first glance the above reasoning may seem to be a little optimistic. It apparently 
overlooks the fact that if we integrate the TDSE over longer and longer times the error on 
the wave function also increases (although it remains bounded because of the unconditional 



stability of the product formula algorithm). In fact it has been shown that in general |17 



||e-*^|V.(0))-^/'"(r)|^(0))||<cA (54) 

where t = mr, suggesting that the loss in accuracy on the wave function may well com- 
pensate for the gain in resolution that we get by using more data in the Fourier transform. 
Fortunately this argument does not apply when we want to determine the eigenvalues as we 
now show. As before we will discuss the real-time algorithm only because the same reasoning 
(but different mathematical proofs) holds for the imaginary-time case. 

Consider the time-step operator (52). Using the fact that any unitary matrix can be 
written as the matrix exponential of a Hermitian matrix we can write 



It is clear that in practice the real-time method yields the spectrum of H{t), not the one 
of H. Therefore the relevant question is: How much do the spectra of H{t) and H differ? 
In Appendix B we give a rigorous proof that the difference between the eigenvalues of H{t) 
and H vanishes as r^. In other words the value of m (or t = mr) has no effect whatsoever 
on the accuracy with which the spectrum can be determined. Therefore the final conclusion 
is that the error in the eigenvalues vanishes as r^/iV where N is the number of data points 
used in the Fourier transform of Tr e~**^(^) . 



IX. RESULTS 

In all our calculations we take J = 1 and h = 0, except for the Ising model in a 
transverse field, where we take h = 0.75. The random numbers Cn,p are generated such that 
the conditions Eqs. (A3) and (A4) are satisfied. We use two different techniques to generate 
these random numbers: 
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1. A uniform random number generator produces {fn,p} and {gn,p} with —1 < /„ p, gn,p < 
1. We then normahze the vector (see Eq. (8)). 

2. The obtained from a two- variable (real and imaginary part) Gaussian random 
number generator and the resulting vector is normalized. 

Both methods satisfy the basic requirements Eqs. (A3) and (A4) but because the first 
samples points out of a 2D-dimensional hypercube and subsequently projects the vector 
onto a sphere, the points are not distributed uniformly over the surface of the unit hyper- 
sphere. The second method is known to generate numbers which are distributed uniformly 
over the hyper-surface. Although the first method does not satisfy all the mathematical 
conditions that lead to the error (14), our numerical experiments with both generators give 
identical results, within statistical errors of course. Also, within the statistical errors, the 
results from the imaginary and real-time algorithm are the same. Therefore we only show 
some representative results as obtained from the real-time algorithm. 

In Fig. |I] we show a typical result for the DOS D{e) of the XY model, the Ising model 
in a transverse field and the mean-field model, all with L = 15 spins and using S = 20 
samples. Because of the very high degeneracy we plotted the DOS for the mean-field model 
on a logarithmic scale. 

In Fig. 0we show the relative error based on Eq. (36) for the three models of various 
size, as obtained from the simulation (symbols). For these figures we used the imaginary- 
time algorithm, because then the statistical error can be related to ca directly (see Eq. (14) 
with A = exp{—i3H)). The theoretical results (lines) for the error estimate, obtained by a 
direct exact numerical evaluation of (14) are shown too. We conclude that for all systems, 
lattice sizes and temperatures there is very good agreement between numerical experiment 
and theory. 

Results for the energy E and specific heat C are presented in Fig. |^ (XY model), ^ (Ising 
model in a transverse field), and ^ (mean-field model). The solid lines represent the exact 
result for the case shown. Simulation data as obtained from S = 5 and S = 20 samples are 
represented by symbols, the estimates of the statistical error by error bars. We see that the 
data are in excellent agreement with the exact results and equally important, the estimate 
for the error captures the deviation from the exact result very well. We also see that in 
general the error decreases with the system size. Both the imaginary and real-time method 
seem to work very well, yielding accurate results for the energy and specific heat of quantum 
spin systems with modest amounts of computational effort. 

X. CONCLUSIONS 

The theoretical analysis presented in this paper gives a solid justification of the remark- 
able efficiency of the real-time equation-of-motion method for computing the distribution of 
all eigenvalues of very large matrices. The real-time method can be used whenever the more 
conventional, Lanczos-like, sparse-matrix techniques can be applied: Memory and CPU re- 
quirements for each iteration (time-step) are roughly the same (depending on the actual 
implementation) for both approaches. 

We do not recommend using the real-time method if one is interested in the smallest 
(or largest) eigenvalue only. Then the Lanczos method is computationally more efficient 
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because it needs less iterations (time-steps) than the real-time approach. However if one 
needs information about all eigenvalues and direct diagonalization is not possible (because 
of memory/CPU-time) there is as yet no alternative to the real-time method. The matrices 
used in this example (up to 32768 x 32768) are not representative in this respect: The 
real-time method has been used to compute the distribution of eigenvalues for matrices of 
dimension 16777216 x 16777216 111. 

Once the eigenvalue distribution is known the thermodynamic quantities directly follow. 
However if one is interested in the accurate determination of the temperature dependence 
of thermodynamic (and static correlation functions) properties but not in the eigenvalue 
distribution itself, the imaginary-time method is by far the most efficient method to compute 
these quantities. For instance the calculation of the thermodynamic properties for f3J = 
0, . . . , 10 of a 15-site spin-1/2 system (i.e. implicitly solving the full 32768 x 32768 eigenvalue 
problem) takes 1410 seconds per sample on a Mobile Pentium HI 500 Mhz system. 

Finally we remark that although we used quantum-spin models to illustrate various 
aspects, there is nothing in the real or imaginary-time method that is specific to the models 
used. The only requirement for these methods to be useful in practice is that the matrix is 
sparse and (very) large. 
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APPENDIX A: EXPECTATION VALUE CALCULATION 

In this appendix we calculate the expectation value of the error squared, as defined in 
Section |T[ By definition we have 



EiYliRA?'] = E 



1 

52 



^ S D 2 

p=l m,n=l 
S D 

p^p'=l k,l,m,n=l 



—D 6m,n E{c^ p,C*i p,) + E{cl^^ p C„ p p, C;*p,) j Al^iAm,n, (Al) 

where p and p' label the realization of the random numbers Cn,p = fn,p + Wn,p- 

First we assume that different realizations p ^ p' are independent implying that 

E{(^*m,p (^n,p Ck,p' Clpi)pj^p' = E{c^ p Cn,p)E{Ck,p' C^^p,). (A2) 

Second we assume that the random numbers are drawn from a probability distribution that 
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is an even function of each variable 

P{fi,pj 9i,pj f2,p, g2,p, ■ ■ ■ , fk,p, gk,p, ■ ■ ■ 1 fD,p, gD,p) 

— P{fl,p, 9l,p, f2,p, g2,p, ■ ■ ■ , —fk,p, 9k,p, ■ ■ ■ , fD,p, gD,p) 

= P{fi,p, 9i,p, f2,p, g2,p, ■ ■ ■ , fk,p, -gk,p, ■ ■ ■ , fD,p, gD,p), (A3) 

and symmetric under interchange of any two variables 

P{fi ,pi ■ ■ ■ ) fi,p) 9i,p) ■ ■ ■ 1 fj,p) 9j,pi • • • 5 fD,pi 9d,p) 
^ P^fl,pi 9l,pi • • • 1 fj,pi 9i,pi • • • 1 fi,pi 9j,pi ■ ■ ■ 1 fD,pi 9d,p) 

,p} 9l,p) • • • 5 9i,pi fi,pi ■ ■ ■ 1 fj,pi 9j,pi ■ ■ ■ 1 fD,pi 9d,p), (A4) 

for all i, j, k = 1, . . . , D. This is most easily done by drawing individual numbers from the 
same even probability distribution i.e. 

D 

P{fl,pi 9l,pi ■ ■ ■ ■> fj,pi gi,pi ■ ■ ■ ■> fi,pi gj,pi ■ ■ ■ ■> fD,pi gD,p) — W P{fn,p)P{gn,p)i 

n,m=l 

where P{x) = P{—x). Normahzing the vector {fi,p, Qi^p, ■ ■ ■ , fD,p, 9d,p) such that 
J2iLi |cn,pP = 1 (for P = 1, • • • , S) does not affect the basic requirements (A3) and (A4). 
Making use of the above properties of P{fi,gi, . . . , fo, go) we find that 

^(Cm,pC„,p) = 5„,^aE{\c„,,pf) = 5„,,nEi\c\'^), (A6) 

where in the last equality we omitted the subscripts of Cm,p to indicate that the expectation 
value does not depend on m or p. An expectation value of a product of two c*'s and two c's 
can be written as 

E{c*^,pC^,pCk,p'Clp,) =(1 - 5py)(5„,„4,/^(|c™,pn^(|c„,,yH 

+ Sp,p'Sm,nSk,l{^ - (^rnk)E{c*„^ 4) 
+ Sp^p'Sm,kSn,li^ - Srn,n)E{c*^, C„ C„ C* ) 
+ 5py6rn,lSnMi^ " 5m,n)-E(c*„ C,„ C„ C^) 
~l~ ^p;p'^'m,l^n,k^'m,nE(Cjyi '-'m) 
= (1 - 5p,p')Sm,n5k,lE{\c\'^)'^ 
+ ^p,p'Sm,nSk,l{^ — Sm,k)E{\ 

Cm,p\ \Ck,p\ ) 

~l~ ^p,p'^m,k^n,l(^^ ^m,n) E (^\c^^p\ \Cn^p\ ) 

+ 3p,p'Sm,l^n,ki^ ~ ^m,n)E{c^ p C^^^ C^ p (^m,p) 

+ 8p,p'5m,l5n,k^m,nE{\Cm,p\'^) ■ (A7) 

Furthermore for m 7^ n wc have 

E{ ,p9m,p 9rn,p )U'n,p + 2^fn,p9n,p- 9n,p)) 

E{fm,pfn,p) ^" '2^E{f„ipfn,p9n,p 

2i£/(/m p(y'„j p/|j p) + 4:E(^f„ipfn.pgm,pgn,p) ~l~ '2^E{fYn,pgm,pgn,p) 
E{9m,pfn,p) '2^E{g„ipfn,pgn,p 

) + E{g 

9n,p) 

=-^(/m,p/n,p) - E{gl^pflp) - E{f^pglp) + E{gl^pglp) 

=0. (A8) 
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By symmetry E{\cmp\'^ Icnpi^) does not depend on m, n or p and the same holds for 



^m,p I 



The fact that the vector (ci^p, . . . , CD^p) is normahzed yields the identities 

E f E \Cn,A = E E{\cnJ') = DE{\c\') = E{1) = 1, 
\n=l / n=l 



(A9) 



and 



E((T.\Cn,A ] = E ^(|Cn,pr|c„,pP) 
y \n=l / J m,n=l 



D D 

= E^(|Cn,p|')+ E (l-5n^,n)^(|c„^C„|2) 
n=l m,n=l 

= DE{\c\^) + D{D - l)E{\c\''\cf) = E{1) = 1, (AlO) 
where c and c' refer to two different complex random variables. Therefore we have 

E{\c\') = 1/D, (All) 



and 



,2,„2> l-DBdcl") 



Substitution into (A7) yields 



E{c*m^pCn,pCk,p' Clp,) =(1 - dpy)Sm,nSk,lD ^ 

-DE{\c\' 
D{D-1) 



+ ^p,p' — {^m,nSk,l{^ — ^m,k) + ^m,k^n,l{^ ~ ^m,n)) 



+ 5p^p'5m,l5n,k5m,nE{\cf). (Al3) 

and the final result for the variance reads 

V ^ ^ S\ D-1 D-1 ' ' 



D-1 



An expression for the fourth moment £'(|c|^) cannot be derived from general properties of 
the probability distribution or normalization of random vector. We can only make progress 
by specifying the former explicitly. As an example we take a probability distribution such 
that for each realization p the random numbers /„ p and (?„ p arc distributed uniformly over 
the surface of a 2D-dimensional sphere of radius 1. This probability distribution can be 
written as 

P{fi,9i,f2,92,... ,fD,9D)^S{f^ + gl + f^ + gl + ... + f^ + gl-l), 

(A15) 
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where we omitted the subscript p because it is irrelevant for what follows. The even moments 

of \cn\ — {fn + QnY^'^ are defined by 



E(|c|2^) 



2M^ _ I-ooifi + gD'^^Ul + 9l + p2+gl + --- + pD + 9l- midgidhdg^ ■ ■ ■ dfndgn 
I-oo S{f! + 9l + --- + pD + 9h- mid9i . . . dfodgo 



(A16) 

It is expedient to introduce an auxiliary integration variable X by 

2M^ I-oo X'^Wi + 9l - X)S{fi + gl + ... + f^ + gl-{l- X))dXdf,dg,df,dg2 . . . dfndgn 



£;(|c|'^) 



!-oo KPi +9! + --- + pD + 9h- midgi . . . dfndgn (^17) 



We can perform the integration over X last and regard (A17) as the M-th moment of the 
variable X with respect to the probability distribution 

I-oo + 9l X)5U1 + gl + ... + fh + gl-{l- X))df\dg,df^dg^ . . . dfndgn 



^^^^ I-oo HP + 9l + --- + pD + 9h- midgi . . . dfndgn 

The calculation of P{X) amounts to computing integrals of the form 



(A18) 



I^{X) = j 5\Y,xl- Xjdxidx2...dxN. (A19) 



Changing to spherical coordinates we have 



yielding 

h{X)hn-2{l - X) 

"^^""^ 

^{D-l){l-X)^-^e{X)d{l-X). (A21) 

The moments E{\c\^^) are given by 

/oo 
X'^P{X)dX 
-00 

= (L> - 1) X^{1 - Xf-^dX 

_ rp)r(i + M) 

" r{D + M) ' ^^^^^ 

and the values of interest to us are 

E(|cr) = l, E{\cr) = ^, mc\') = ^J^^y (A23) 

where the first two results provide some check on the procedure used. Substituting (A23) 
into (A14) yields 

E{\TrRAf)= ^^^^[^ I . (A24) 
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APPENDIX B: ERROR BOUNDS 



Here we prove that the difference between the eigenvalues of the Hermitian matrix A + B 
and those obtained from the approximate time-evolution exp (2^4/2) exp(zi?) exp(2A/2) [z = 
—IT, — r) is bounded by r^. In the following we assume A and B are Hermitian matrices 
and take r a real, non-negative number. We start with the imaginary-time case. 

We define the difference -R(t) by 

R{t) =e^(^+^) - e^^/2grBgrA/2 

=\ r dX l^dfi rrfi/e^^/'e^^{e-"^[25,[A,5]]( 
4 JO JO JO 



,uB 



a well-known result |21|. We have ||23[ 

mr)\\<l 



(Bl) 



< 



rdX I dfi rrfi/e^^/2g("-'^)^[2i?,[A,S]]e^«e^^/2e(^-^)(^+^) 

JO JO JO 

^ dX f\fi rdz/e^^/2gABg.A/2M^ [A,B]]e 

JO JO 

\ r dX f'dfi rtiz/e^ll^ll/V^-^)ll^ll||[25,[A,5]]||e^ll^lle 
4 Jo JO JO 



1 



{X-u)A/2^{t-X){A+B) 



A||A||/2g(r-A)(||A|| + ||B||) 



■4 

+ i rfA t^/i f c^^^e^ll^ll/ Vll^lle'^ll^ll/^l I [A [A S]] | |e(^-)ll^ll/V-^)(ll^ll+ll^ll) 











-rV(ll^ll+ll^ll)(||[A, [A,B]]\\ + \\[2B,[AB]]\\), 



(B2) 



and 



\\Ri-r)\\<- 



+ 
1 

'A 



dX f\fi rdi/e"^/2e("-'^)^[2i?,[A,5]]e'^^e"^/2g(-r-A)(A+B) 



JO 

A rfJ- 



JO 
A 



dX dfx diye^^/^e^''e''^/^[A,[A,B]]e- — '"e 







(A-^)A/2^(-r-A)(A+B) 
-T+A)(A+B) 



1 

+ 4 



JO JO JO 

rdX l^dfi r dve~'^''/\'^^e-'^/^[A, [A,B]]e 

JO JO JO 



{-\+v)A/2 ^{-T+\){A+B) 



<i fdX f'df, rdue'\\''\\/'e^'-^^\\''%2B,[A,B]]\^^^^^^^ 
4 JO JO JO 

+ - TdX I'dfi rdue'\\''\\/'e'\\V''\\/'\\[A,[A,^ 
4 Jo Jo Jo 

^ ' " " " " (B3) 



.3^r(||A|| + ||i?||) 
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r-e ^-{\\[A, [A,B]]\\ + \\[2B,[AB]]\\). 



Hence the bound in -R(r) does not depend on the sign of r so that we can write 

\\Rir)\\ <s|r|V^I(ll^ll+ll^ll), (B4) 



17 



where 

s^^\\[A[AB]]\\ + \\[2B,[A,B]]\\. (B5) 

For real r we have 

^rA/2^rB^rA/2^^rC(r)^ (B6) 

where C(r) is Hermitian. Clearly we have 

e-(^+fi) _e-^M =i?(r). (B7) 

We already have an upperbound on R{t) and now want to use this knowledge to put an 
upperbound on the difference in eigenvalues of C(r) and A+B. In general, for two Hermitian 
matrices U and V with eigenvalues {«„} and respectively, both sets sorted in non- 
decreasing order, we have p| 

|Mn < ■^ll, Vn. (B8) 

Denoting the eigenvalues of A + B and C(r) by Xn^O) and Xn(r) respectively, combining 
Eq. (B4) and (B8) yields 

|grx„(o) _ g™„(r)| ^ slrl^el^K'l^ll+'l^'l). (B9) 

To find an upperbound on |x„(0) — a;„(r)| we first assume that x„(0) < x„(r) and take 
r > 0. It follows from (B9) that 

For X > 0, — 1 > X and we have — a;„(0) < | |y4 + i?| | < | |y4| | + | |. Hence we find 

Xn{r) - Xn{0) < sr^e^^^ll^ll+ll^ll). (Bll) 

An upperbound on the difference in the eigenvalues between C (r) and A + B can equally 
well be derived by considering the inverse of the exact and approximate time-evolution 
operator (B6). This is useful for the case x„(0) > x„(r): Instead of using (B7) we start 
from exp(— r(A + B)) — exp(— rC(— r)) = R{—t) (r > 0). Note that the set of eigenvalues 
of a matrix and its inverse are the same and that the matrices we are considering here, i.e. 
matrix exponentials, are nonsingular. Making use of Eq. (B4) for R{—t) gives 



e 



-rx„(0) _ g-™„(r)| ^ _g|^|3g|r|(||A|| + ||B||)^ |-gj_2) 



and proceeding as before we find 

r(x„(0) - x„(r)) < e^(-"(o)-xn(r)) _ ^ < 3^3g2r(||A||+||B||)^ (3^3) 

Putting the two cases together we finally have 

|x„(r)-a;„(0)|<srV^(ll^ll+ll^ll). (B14) 
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Clearly (B14) proves that the differences in the eigenvalues of A + B and C(r) vanish as r^. 

We now consider the case of the real-time algorithm {z = —ir). For Hermitian matrices 
A and B the matrix exponentials are unitary matrices and hence their norm equals one. 



This simplifies the derivation of the upperbound on R^—ir). One finds 

\\R{-^r)\\E < s\t\\ (B15) 

where = TrA'''^! denotes the Euclidean norm of the matrix A [0. In general the 

eigenvalues of a unitary matrix are complex valued and therefore the strategy adopted 
above to use the bound on R{t) to set a bound on the difference of the eigenvalues no longer 
works. Instead we invoke the Wielandt-Hoffman theorem p3: 



If U and V are normal matrices with eigenvalues ui and Vi respectively, then there exists a 
suitable rearrangement (a permutation g of the numbers 1, . . . ,n) of the eigenvalues so that 

J:\u,-v,^,)\'<\\U-V\\1. (B16) 

Let U and V denote the exact and approximate real-time evolution operators respectively. 
The eigenvalues of A + B and C(r) are x„(0) and x„(r) respectively. All the XnS and a;„(r)'s 
are real numbers. According to the Wielandt-Hoffman theorem 

N 

i=i 

where |/j(r) = a;p(j)(r), g being the permutation such that inequality (B17) is satisfied. We 
see that Eq. (B17) only depends on (rXj(O) mod 27i) and {ryjij) mod 27r), but so does 
the DOS (see Eq. (|16D). Since the inequality (B17) and the DOS only depend on these 
"angles" modulo 27r, there is no loss of generality if we make the choice 

0< |r(a;,(0)-y,(r))| <7r. (B18) 

Rewriting the sum in (B17) we have 

N N 

^ |e^--i(o) - e'"?^iM|2 ^ - 2 cos(r (x,-(0) - %(r)))) 
i=i i=i 

N 

= 45:sin2(r/2(x,(0)-i/,(r))). (B19) 

As we have 

sin^ X < for < |x| < 7r/2, (B20) 

the restriction Eq. (B18) allows us to write 

N 2 2 

$:(x,(0)-y,(r))^<^r^ (B21) 
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implying 



\xA^)-yAr)\<-T\ (B22) 

In summary, we have shown that in the real-time case there exists a permutation of the 
approximate eigenvalues such that the difference with the exact ones vanishes as r^. 

Finally we note that both upperbounds (B22) and (B14) hold for arbitrary Hermitian 
matrices A and B and are therefore rather weak. Except for the fact that they provide 
a sound theoretical justification for the real and imaginary-time method, they are of little 
practical value. 
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FIG. 1. The density of states (DOS) as obtained from the real-time algorithm for spin chains 

of length L = 15 and for S* = 20 random initial states. Left: XY model; middle: Ising model in a 
transverse field; right: Mean- field model. For the mean- field model a logarithmic scale was used 
to show the highly-degenerate spectrum more clearly. 





FIG. 2. The relative error 5Z/Z (see Eq. (36)) on a logarithmic scale as a function of temper- 
ature T =1/13 and for various system sizes. Left: XY model; middle: Ising model in a transverse 
field; right: Mean-field model. Solid lines: ea (with A = e~^^ , see Eq. (14)) for L = 6; dashed 
lines: ca for L = 10; dash-dotted line: eA for L = 15. Crosses: Simulation data for S = 20 and 
L = 6; squares: Simulation data for S = 20 and L = 10; circles: Simulation data for S = 20 and 
L= 15. 
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FIG. 3. Energy (top) and specific heat (bottom) of the XY-model (see (39)) with A = 1, 
h = and J = 1. Left: L = 6; middle: L = 10; right: L = 15. Sohd hnes: Exact result, crosses: 
Simulation data using S = 5 samples; squares: Simulation data using S = 20 samples. Error bars: 
One standard deviation. 




FIG. 4. Energy (top) and specific heat (bottom) of the Ising model in a transverse field 
(sec (39)) with A = 0, h = 0.75 and J = 1. Left: L = 6; middle: L = 10; right: L = 15. Solid 
lines: Exact result, crosses: Simulation data using S = 5 samples; squares: Simulation data using 
S = 20 samples. Error bars: One standard deviation. 
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FIG. 5. Energy (top) and specific heat (bottom) of tlie mean-field model (see (45)) witli J = 1 
and h = 0. Left: L = 6; middle: L = 10; right: L = 15. Solid lines: Exact result, crosses: 
Simulation data using S = 5 samples; squares: Simulation data using S = 20 samples. Error bars: 
One standard deviation. 
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